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Modified  rhoJe!'.ky  Decomposition 

LANOUAGE 

ISO  Fortran 

DESCRIPTION  AND  PURPOSE 

let  Y(l),  ...,  Y(T)  be  a  sample  realization  of  a  mixed  autoregressive 
moving  .average  time  series  !Y(t),  t  =  0,  +  1,  ...1  of  order  (p  ,q)  (denoted 
ARM/Mp,q}  ).  Thus 

P  q 

\  »(,1)Y(t-,i)  =  y  B(k)t(t-k)  ,  t  =  0,  +  1,  ... 
j=0  k=0 

for  (on.st.ints  p,q  ,  <i(0)  =  B(0)  “  1  ,  i(l),  ....  ri(p),  B(l),  ...,  S(q)  , 

where  (  (•)  is  a  white  noise  serie.s  of  zero  mean  uncorrelated  random  variables 

2  .P  j 

having  variance  0  .  The  zeros  of  the  complex  polynomials  g(z)  = 

r"'  k 

and  h(z)  =  /j^_Qp(k)z  are  assumed  to  be  outside  the  unit  circle. 

Subroutine  MXPD  calculates  exact,  memory  h,  horizon  t,  minimum  mean 
square  error  linear  predictors  Y(t+h|t)  and  (optionally)  prediction  variances 
I  of  Y(t+h)  given  Y(l).  ....  Y(t)  foi  h  =  h^ ,  ....  hj,  and  t  .... 

Thu;;  Y(t+h|t)  l.s  that  linear  combination  of  Y(l),  ...,  Y(t)  closest  to 
Y(t+h)  in  the  mean  square  sense  and  o^  =  E{Y(t+h)  -Y(t+h|t))  is  the 
attained  minimum  mean  square  prediction  variance. 

If  q  ~  0,  Y  is  a  pure  autoregressive  process  of  order  p  (AR(p)),  while 
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if  p  =  0  ,  Y  is  a  pure  moving  average  process  of  order  q  (MA(q)).  Subroutines 
ARPD  and  MAPD  calculate  Y(t+h|t)  and  (optionally)  for  the,  AR{p)  and 

MA(q)  cases  respectively.  Separate  subroutines  are  used  for  these  c.iscs 
due  to  significant  simplifications  in  their  a  I gor i t hn s  from  the  gencrn! 
ARMA(p,q)  case. 


NUMERICAL  METHOD 

I'ho  method  Is  basi'd  on  numerous  spc'cl.tl  ji  >  ■;  !•  r  t  I  I's  of  the  nioii  1  f  U  ct 

Cholesky  decomposition  (MCD,  see  Wilkinson  (lh67'/)  o."  the  autocovar !  ance 

matrix  of  an  ARMA(p,q)  process.  (See  Newton  and  Pag.ano  (1930)  for  details), 

Let  [Aj^j  denote  the  (i,j)th  element  of  ,c  .n.u  r . A  ind  de.i'ine 

A^  r  T0EPL{a(0),  ....  a(K-l)}  to  be  the  (K^K)  synine.trit  Toeplitz  matrix 

having  fA^I^j  “  a(|i"j|)  •  het  A^,  =  Lj^Dj^L^  ,  K  >  i  he  the  MCD  of  the 

symmetric  positive  definite  nested  matrices  Aj^ ,  A,,  i._e. 

=  »  I-jf  is  a  unit  lower  t :  .  an  (K'K)  mitrix,  and 

D  is  a  diagonal  matrix.  Then  the  se(|uences  o;  mat  ■'ces  I,  D  ,  , 

K '  K.  \ 

and  are  also  nested  and  we  write  (for  example;  ae  ll,j)tli  element 


of  any  L^^  for  K  >  max  (1,i)  as  ll-ljj  • 


Thus  (D) 


11 


[A] 


11  ’ 


[..I,,  -  IIAI,, 


ir' 

),  1 1'  I  j  1  D I  g  g  f  I  ■  I  ,  ■ 


i-l 


(nil,  =  fAl,, 


f  -  1 


i  •  K 


We  define  the  following  quantities  for  K  1 

a)  r  =  TOEPI,  IR./0),  -  K..(K-I)1  -  1...  .  I  ' 

/• ,  K  /.  A  /  ,  K  /  ,  K  ,  K 


Rj,(v)  -  E|7.(t)Z,(t+v)  ( 


0,  +  1.  . 


wile  re 
nul  Z(-)  .m  AR(p) 


-J- 


process  linvlnR  coefflclentK  rt(l),  <Ti(p)  and  white  noise 

variance  . 


b) 

II 

X  ' 

(X(]).  .... 

XCK))"^  =  -'zIk 

Y  ,  where 

K 

T 

Y  = 
-K 

(Y(l),  .... 

Y(K))  . 

c) 

r 

X,K 

-  - 

T 

4.k'^x,k''x,k  ■ 

Note  that  I'  „  = 

A  y  K 

d) 

II 

0)^ 

(e(l),  .... 

e(K))^  ■  I-xIk 

4  • 

-1 


,-T 


Then 


Y(t+h|t)  =  X(t+h|t)  -  I  oi(j)Y(t+h-jit) 

j  =  l 

2  2 

^t,h  "  |^^g^^Z^X^t+h,t+h-k^’^X^t+h-k,t+h-k 


Jo 


where  Y(t+h-j|t)  =  Y(t+h-j)  If  j  >  h,  and 
X(t+h|t)  = 


,t+h-k|^^^^.^J^_^^^4^t+h,^^4^^,t+h-kj  I 


q 

X  ^ . . 

k=h 


,  h  >  q 


Thus 


2 

to  obtain  Y(t+h|t)  and  ^  for  h  =  h^^,  ...,  h2  and  t  =  tj^,  ....  t^ 


,-l 


4.t2+h2  *  4;t2+h2  ’  4,t2+h2  4.t2+h2  '  Significant  reductions 

In  computing  time  and  storage  requirements  in  obtaining  the  elements  of 
these  matrices  are  afforded  by  noting: 


NOTE _J.  Comput ing  L, 


-1 


2^**2 


The  row  of  I-,',,  is  given  by 


'I- 


(1,  , 


....  Oj_j(l),  1,  Q^.j)  .  2  £  j  _<  p 


(0 


]_p_l  •  n(p),  «(1),  1,  p+1  £  j  1  K 


(1) 


(2) 


where  “  «(^)  If  k  >  p  and 

a  .  (1 )  "  -  - ^  ^  - - - 

i-<’‘j+i(j+i) 


;  i  -  1 . J  p 


llius  there  are  only  p(p+l)/2  distinct  nonzero, nonone  elements  In  I. 
K  >  p  +  1  and  X(l)  =  Y(l)  while 


-1 

Z.K  ’ 


X(j)  = 


Y(j)  +  I  ,(i)Y(,i-£)  ,  j 
f=l 


Y(j)  +  I  a(^)Y(j-^) 
£=1 


=  2,  ....  p 


(3) 


j  >  P 


NOTE 


2  Computing 

•  L  ^ 


Let  Y(0)  =  1,  y(1),  y(2),  ...  be  the  coefficients  of  the  MA(«>) 
representation  of  Z(*)  ,  1 .e . 

j-1 

Y(j)  =  -  1  a(j-i)y(i)  ,  j  >  1  ‘ 

£=max(0, j-p) 


Then 


^4^p+j,p+j-£  '  Y(i) 


r=l  P+J-r,£ 


,  0  <  ^  <  j  <  K-p 

,  ^  =  1 . P-1,  j  ±  1 


L„  =  (L  ‘  ^  11 

Z ,  p  Z ,  p 


u  ^4^1r^4^r,J-k’ 

r=1-V-»-l 


P-1 


lim 


1  im 

K - #  00 


J  "'zIk.j  ■» 


K 


"-A.. 


^  =  R^(0)/o^ 


(4) 

(5) 

(6) 

k<J<p 

(7) 

(8) 

(9) 


—  )— 


y  Y  (k)  =  R_(0)/a^ 
k=0  ^ 


i/O)  =  a!  II  (l-rt,(j)) 

'i  =  l  J 


By  (5)  <nnd  (6)  the  distinct  elements  of  ,,  are  contained  in  Its 

Z,t-+h- 

T  ^  ^  T 

first  p  columns,  the  last  of  which  Is  (0^  1,  y(1)»  •••»  Y  ^  ^ 

Thus  MXPD  finds  the  elements  of  L„  ,,  by: 

Z,t2+h2 

i)  finding  L„  =  (lZ^  )”^  ,  11)  calculating  Y(l) . Y(M,-p) 

7, ,  p  Z ,  p  1 

(by  (4))  where  Mj^-p  is  the  smallest  integer  such  that 


y  Y  (j)  -  R„(0)/o^  <  6  ,  (12) 

j=0 

ill)  use  (6)  to  obtain  remaining  elements  of  rows  p+1,  of 

first  p-1  columns  of  L  .  .  Note  that  (10)  says  that  such  an  M. 

exists  (which  hopefully  is  much  less  than  t2+h2)  while  (8)  and  (9)  say 

that  all  further  elements  of  the  first  p  columns  of  L„  .  are  arbitrarily 

Z,t2+h2 

close  to  zero.  Thus  there  are  M^^p  elements  to  calculate  and  store  for 

4,t2+h2  ■ 


NOTE  3 


Computing  L 


X,t2+h2  ,  X,t2+h2 


To  compute  _  ..  and  „  ..  we  note  that  (defining  a_(0)  =  1) 
X,t2+h2  A,t2+h2  u 

[r„]..  =  r  1  n4_i(J-ni)  I  a  Ai-DRyU-m),  IJ^l 

^  I  m=max(l,J-P)  ^  ^=max(l,l-p) 


.  i.j  >  p.  |i-j  I  1  q 


,  1  <  J  £  p,  1  >  p,  and 
1  -  j  >  q  or  if 


-f)- 


RyCv)  =  I  e(k)B(k+v)  ,  v  =  0,  . . . ,  q  .  (1^) 

^  k=0 

Thus  the  distinct  elements  of  ^  are  Rj^CO)  >  •  ■  • » 
elements  In  the  first  p  +  q  rows  and  p  columns.  Further  only  Ry(0), 

R^Cp+q)  are  required  to  obtain  .  These  elements  are  obtained  via 

subroutine  MXCV  by  solving  first  for  v  =  0,  max(p,q) 

^  a(j)RY(v-J)  =  J  e(k)R^  (v-k)  =  0  ,  V  >  q  (15) 

j=0  -  k=v 

'y 

where  Ry£(v)  -  E{ Y (t) c (t+v) }  =  if  v  >  0  ,  where  ^  Is  the  Kronecker  delta. 


2  mln(v,p) 

R^^C-v)  =  e(v)a  -  I  a(j)RY^(-v+j)  ,  v  =  1,  q  .  (16) 

j=l 

Then  RY(max(p,q)+l) ,  RY(p+q)  are  obtained  by  (15)  for  v  =  max(p,q)+l, 

....  p+q  . 

To  obtain  ^  ..  and  D  we  note  that  the  pattern  of  zeros  In 

A  j  c  2  ’  n  2  “ ♦ ^  2^^  2 

L„  .  Is  the  same  as  that  In  the  lower  triangle  of  r„  ,,  and  that  th« 
X,t2+h2  X,t2+h2 

required  elements  of  F  can  be  calculated  as  needed  (by  (13))  without 

X,  p+q 

storing  them. 

Thus  L„  Is  calculated  and  stored  In  one  matrix.  To  obtain  the  q 

X,p+q 

nonzero  nonone,  elements  of  the  rows  of  the  rest  of  L„  ,  we  have: 

A  j  1 2*' ^2 


lln  141,^,.,  -  li(J) 


*  j  ~  1*  •••»  ^ 


lim  [DxV.K 


Thus  rows  p  +  q  +  1,  ...  are  calculated  and  stored  In  a  matrix  having  q 


columns  until  the  elements  of  I.  ,  U  have  converged  (at  row  M_  say)  ,  i.e. 

XX  ^ 

I  I4l„  I  -  «(q-J«)l  '  ».  II4I„  „  -  "'1  '  J  -  V’ . “2-‘ 

2 ,  2  2 


-H- 


iv)  Find  and  M  +p^ij  ’  *  *  1  1  )  1  P  by  (7),  (6), 

and  (12)  (stored  in  the  integer  MONK  and  (Mj^xp)  matrix  ALZ)  . 

v)  Find  R  (0),  RY(q)  by  (14)  via  subroutine  MACV  (stored  in 

A  A 

constant  RXO  and  q-vector  RX) . 

vi)  Find  L„  .  and  1)„  .  (stored  in  the  (p+q)x(p+q)  matrix  ALXl 

X , p+q  X , p+q 

and  in  the  first  (p+q)  elements  of  the  M^-vector  DX) . 

vil)  Find  and  ^  ^  P+q+1 .  •••»  ^2, 

j  =  i-q,  ....  i- 1  (stored  in  the  integer  MTWO  and  the  (M2Xq) 

matrix  ALX2  and  the  rest  of  the  (M2)-vector  DX) . 

viii)  Find  X  ,  e  by  (3)  and  (18)  (stored  in  the  t.-vectoisX  and  E)  . 
'^2  '^2 

lx)  Find  Y(t+h|t)  ,  for  h  =  hj^,  •••>  h2  ,  t  =  t^,  t2  (stored 

in  the  (t2-t  j^+1)  (h2-hj^+l)-vectors  YPD  where  Y(t+h|t)  =  YPD((t-tj^) 

(h2-h^+l)+(h-h^+l))). 

2 

x)  Find  (optionally)  j.  which  Is  stored  like  YPD  in  a  (t2-tj^+l) 
(h2-hi+l)-vector  PVAR  . 


NOTE  6  Taking  advantage  of  convergence 

To  take  advantage  of  the  convergence  in  and  the  user 

specifies  an  absolute  convergence  criterion  6  and  Integers  IROWSl, 

1ROWS2  as  the  row  DIMENSION  of  arrays  ALZ  and  AiL.X2 ,  DX  respectively. 

Thus  IROWSl  and  IROWS2  must  be  chosen  to  exceed  what  can  reasonably  be 
expected  to  be  and  M2  respectively  or  else  a  nonconvergence  failure 
Indiiator  will  be  returned  in  IFADI.T.  Of  course  if  IROWSl  or  1ROWS2  are 
given  the  value  i2^^'2  ^''cn  convergence  need  not  be  reached  for  the  algorithm 
to  finis  properly.  Setting  IROWSl  or  IR()WS2  smaller  than  t2^^2  *bc 

possibility  of  obtaining  predictors  for  long  time  series  with  a  minimum 
amount  of  required  storage. 


-9- 


Further,  if  the  ^  are  not  to  be  calctilateci ,  the  nuatrix  ALZ  Is  not 
needed  and  I  ROWS  1  can  be  set  equal  to  1. 

NOTE  7  Algorithm  for  ARPD 

For  q  =  0  we  have  F,  .  =  ,  r  =  D  L  =  I  D 

Z,K  Y,K  'x.K  "^Z.K  ^X,K  ’  ‘^X.K 

=  .  Thus 


Y(t+h|t)  =  -  ^  a(j)Y(t+h-j |t) 

j  =  l 


2  2 
=  I  IL  ] 


t,h 


[l)J 


k=0 


Z't+h,t+h-k  ^  Z't+h-k,t+h-k 


2  9 

=  0  [  Y  (k)  if  t  >  p 

k=0 


NOTE  8  Algorithm  for  MAPD 


TOEPI 


For  p  -  0  ve  hove  -  1,^  e„d  . 

fR^(O),  R^(q),  0,  ....  0}  ,  (so  that  calculation  of  ALXl  ip 


avoided),  K  -K  '  -K  ' 


Y(t+h|t)  = 


X  ^’^X^t+h.t+h-k  ^<t+h-k) 


t,h 


B  (k)e(t+h-k) 

k=h 

0 

h-1 


I  IL 


k=0 


X ' t+h , t+h-k  ^ ”x^  t+h-k , t+h-k 


9  1>-1  a 

I  (k) 

k=0 


,  t+h  <  M,, 


,  t+h  >  M,, 


,  h  >  q 


,  t  <  M„ 


,  t  >  M, 
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STRUCTURK 

SUBROUTINE  MXPD  (NP,  NQ,  AI.PUA,  BETA,  SICSQ,  Y,  lOPT,  NTl ,  NT2 ,  NHl,  NH2 ,  NYPD, 
NPVAR,  NPPNH2,  IROWSl,  IR0WS2,  IROWS  1,  DEI.,  RYE,  RYEO,  RY ,  RYO ,  IP,  YWK,  RZO, 
RX,  RXO,  ALZl,  ALZ,  ALXl,  ALX2 ,  DX,  MONE,  MTWO,  X,  E,  YPD,  PVAR,  IFAUl.T) . 

Forma  I  parameter's 


NP 

Integer 

input : 

order  of  AR  part  of  model 

NQ 

Integer 

input : 

order  of  MA  part  of  model 

ALPHA 

Real  Array  (NP) 

i nput : 

coefficients  of  AR  part  of  model 

BETA 

Real  Array  (NQ) 

i nput : 

coefficients  of  MA  part  of  model  • 

SIGSQ 

Real 

input : 

variance  of  white  noise  in  model 

>• 

Real  Array  (NT2) 

input ; 

data  vector 

lOPT 

Integer 

input : 
1 

0 

option  switch  equal  to: 

if  both  predictors  and  variances 

to  be  calculated 

if  only  predictors  desired 

NTl 

Integer 

i uput ; 

tj  (first  memory) 

NT2 

Integer 

input : 

t^  (last  memory) 

NHl 

Integer 

i nput : 

h^  (first  horizon) 

NH2 

Integer 

input ; 

h,^  (last  horizon) 

NYPD 

Integer. 

input : 

(NT2-NT1+1) (NH2-NH1+1) 

NPVAR 

Integer 

input : 

same  as  NYPD  if  lOPT  =  1,  >1  if  IOPT=0 

NPPNH2 

Integer 

Input : 

NP  +  NH2 

IROWSl 

Integer 

input : 

row  dimension  of  ALZ  in  calling  pro¬ 
gram  (^NP+Z  if  lOPT  =1,  ^1  if  IOPT=0) 
see  note  6  above 

I  ROWS 2 

Integer 

Input : 

row  dimension  of  ALX2,  DX  in  calling 
program  (>NP+NQ-tl)  See  note  6  above. 

I  ROWS  3 

Integer 

input : 

row  dimension  of  RY,  IP,  ALZI,  ALXl 
in  calling  program  (>NP+NQ) 

DEL 

Real 

input : 

absolute  convergence  criterion  (see 
(12)  and  (17)). 

RYE 

Rea]  Array  (NQ) 

output : 

Ry,(-1) . 

RVEO 

Real 

output : 

R^,(0) 

RY 

Real  Array  (NP+NQ) 

output ; 

. Ry(P+q) 

HYO 

Real 

output : 

Ry(0) 

I  !• 

Integer  Array  (IROWSl) 

work : 

VWK 

Real  Array  (NPPNII2) 

work : 

IJ/.O 

Rea  1 

out  put  : 

Ry(n) 

1  X 

Real  Array  (NQ) 

out  put  ; 

Ky(l) . Rx(q) 

HXO 

Real 

out  put : 

Ry(0) 

\L7. 1 

Real  Array (IROWSl,! ROWSJ) 

out  put : 

L'.p 

ALZ 

Real  Array  (IROWSl,  NP) 

output : 

if  lOPT  =  1,  ALZ  contains 

i  =  1 . p.  1  =  I . , 

iC  IOP7'  =  0,  AI,/  Is  not  used. 
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AI.X  1 

AI.X2 

Kc.-il  Army 

Rcvil  Array 

(  IROWP.'l,  IROWS  J) 

(  IR()WS2,NQ) 

out  |>iil  : 

output : 

''X,p  tq 

|l,j^  1  ^  .  ,  1  =  p+q+1,  ,  .  .  ,  , 

liX 

Kc-al  Array 

(1  ROWS 2) 

out  put : 

i  ^  l-'l . 1-1- 

(n^l,  .  i  =  1,  M2 

MONK 

hUcRor 

out  pul  : 

Mj  (see  note  6  above) 

M  TWO 

1  nt  <■(>(.■  r 

output : 

(see  note  6  above) 

X 

Real  Array 

(NT2)  • 

output : 

Xd) . X(t2) 

1' 

Real  Array 

(NT2) 

output : 

e(l) . e(t2) 

YI’O 

Real  Array 

(NYPD) 

output : 

((Y(t+htt),  h  =  h^,  ....  h2)  . 

- 

t  =  tj^,  . .  . ,  t2) 

PVAR 

Real  Array 

(NPVAR) 

output : 

if  lOPT  =  1,  ((0^  h’  ° 

if  lOPT  =  0,  PVAR  not  used. 


1  FAULT  Integer 

output:  failure  indicator 

r  v'luv''  in.iicationp 

Value  of  1  FAULT 

Meaning 

0, 

no  failure 

1 

NP<1,  NQ^l,  or  lOPT  not  0  or  1 

2 

NTlrNP+NQ  or  NT1>NT2 

3 

NH1<1  or  NH1>NH2 

k 

NYPD<(NT2-NT1+1) (NH2-NH1+1)  or  NPVAR<1 
or  lOPT  =  1  and  NPVAR< (NT2-NT1+1) 
(N112-NH1+1) 

5 

IROWSl-NP+2  and  lOPT  =  1  or  IR0WS1<1 
or  NPPN1I2- NP+NU2 

6 

lROWS2<NP4NQ+l 

7 

IROWSl-'NP+NO 

8 

SK.SIL  0 

9 

Singular  matrix  in  iubroutine  MXCV 

10 

An  0  1  ( j  )  1  (see (2)  ) 

11 

lOPT  =  1.  1R0WR1'NT2+NH2-NP,  and 
convergence  not  reached. 

12 

Non post  live  (0^1,,  encountered 

A  11 

13 

I ROWS2' NT2+NH2-NP-NQ  and  convergence 
not  reached . 

SllBROUriNE  ARPD  (NP, 

ALPHA,  SIGSQ, 

Y,  lOPT, 

Nil,  NT2,  NHl ,  NH2,  NYPD,  NPPNH2 

YWK,  GAM,  YPD,  PVAR, 

I  FAULT) . 

/■'(irma! 

pavnmetnrs 

NP 

Integer 

input 

order  of  AR  model 

ALPHA 

Real  Array 

(NP) 

i  nput 

coefficients  of  AR  model 

SiGSQ 

Real 

i  nput 

variance  of  white  noise 

Y 

Real  Array 

(NT2) 

input 

data  vector 
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lOPT 

Integer 

Input ; 

option  swltcli  equal  to 

1 

if  both  predictors  and  variances 

to  be  calculated 

0 

If  only  predictors  desired 

NTl 

Integer 

Input : 

tj  (first  memory) 

NT2 

Integer 

input : 

(last  memory) 

NUl 

Integer 

input : 

(first  horizon) 

NI12 

Integer 

input : 

h^  (last  horizon) 

NYPD 

Integer 

input : 

(NT2-NT1+1) (NH2-NH1+1) 

NPPNH2 

Integer 

input : 

NP  +  NH2 

YWK 

Real  Array  (NPPN112) 

workspace ; 

r.AM 

Real  Array  (MPPNH2) 

out  put : 

y(1) . y(NP+NH2) 

VPD 

Real  Array  (NYPD) 

output : 

((Y(t+h| t).  h  =  h^,  ....  h2). 

PVAR 

Real  Array  (NH2) 

output : 

if  lOPT  =1,  o^h=l,  ...,h 
t ,  n 

if  lOPT  =  0,  PVAR  not  used. 

I FAULT 

Integer 

output : 

failure  indicator 

Failure  Indiaations 

Value  of  I FAULT 

Meaning 

0 

no  failure 

1 

NP  <  1 

2 

NTl  <  NP  or  NTl  >  NT2 

3 

NUl  <  1  or  NHl  >  NH2 

4 

SIGSQ  ^  0 

5 

lOPT  not  0  or  1  or  NYPD 

<  (NT2-NT1+1) 

(NH2-NH1+1)  or  NPPNH2 

<  NP+NH2 

SUBROUT  INF.  MAPD  (NQ, 

BETA,  SIGSQ,  Y 

,  lOPT, 

NTl,  NT2,  NHl,  NH2, 

NYPD,  NPVAR, 

IRONS, 

DF.L,  RX,  RXO, 

DX,  ALX,  MTWO, 

E,  YPD, 

PVAR,  I FAULT) 

'-'ormaL  parameters 

NQ 

Integer 

input ; 

order  of  MA  model 

BF.TA 

Real  Array 

(NQ) 

input ; 

coefficients  of  MA 

model 

SIGSQ 

Real 

input : 

variance  of  white  noise 

Y 

Real  Array 

(NT2) 

Input : 

data  vector 

[OPT 

Integer 

Input : 

option  switch  equal 

to: 

1 

if  both  predictors 

and  variances 

to  be  calculated 

0 

If  only  predictors 

desired 

NTl 

Integer 

input : 

t^  (first  memory) 

NT  2 

Integer 

input : 

t2  (last  memory) 

NUl 

Integer 

Input : 

hj  (f Irst  hoi .zon) 

NH2 

Integer 

input : 

h2  (last  horizon) 

NHl 

NH2 
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NVPO 

Integer 

input  : 

C.'T.'-NI  Mil  (NII2-NI11+1) 

Nl’VAR 

Integer 

i nput : 

Same  as  NYPD  it  lOPT  =  1, 

•  1  if  I  OPT  -  0. 

I  Rows 

I ntegcr 

Input : 

row  dimension  of  AI.X  and  I)X  in 
i.illing  program  (see  notes  6  and 

H  alioviO 

OKh 

Ron  1 

i nput  : 

alisolute  convergence  criterion 

RX 

Real  Array 

(NQ) 

out  put  : 

R.,(I) . R^(q) 

RXO 

Real 

out  pul  : 

R,.  (0) 

OX 

Rea]  Array 

(IRONS) 

out  put  : 

1  =  1 . M2 

A1,X 

Real  Arrav 

( IRONS, NQ) 

out  put  : 

,,  i  =  1,  ....  M2, 
i  =  i-q,  . . .  ,  1-1 

MI'NO 

Integer 

out  put : 

M  ,  (see  notes  6  and  8  above) 

E 

Real  Arrav 

(NT2) 

out  put : 

e  ( 1 ) . e  ( 1 2 ) 

Yl’D 

Real  Array 

(NYPD) 

output : 

((Y(t+hlt),  h  =  h^ . h2), 

t  -  t^,  12) 

PVAR 

Real  Array 

(NPVAR) 

output : 

if  lOPT  =  1,  ((0^  h  =  hj^ . 

112),  t  =  t^,  12),  if  lOPT  =  0 

PVAR  not  used 

J  FAllLT 

Integer 

output ; 

failure  indicator 

F'sHure  Indiaationr, 

Value  of  IFAULT  Meaning 


0 

1 

2 
3 
U 


5 

6 

7 

8 


no  failure 

NQ  •'  1  or  lOPT  not  0  or  1 
NTl  ^  NQ+1  or  NTl  ■>  NT2 
NHl  -  1  or  NHl  >  NH2 
NYPI)  <  (NT2-NT1+1)  (NH2-NH1+1) 
or  NPVAR  <  1  or  lOPT  =  1 
and  NPVAR  <  (NT2-NT1+1) (NH2-NH1t1) 

IRONS  <  2 
STGSQ  <  0 

nonposltivo  [D  )..  encountered 

A  XX 

IRONS  <  MT2+NH2  and  convergence  not  reached 


A'.ixi  linr'j  al.ijorithrir, 

SUBROUTINK  MACV  (NQ,  BKTA,  SIGSQ,  RX,  RXO,  I FAUhT) 

; '  -rmiiL  pa  rame  t  c  rr, 


NO 

I’.l'TA 

si(;sQ 


r lit  eger 

Rc.il  Array  (N(^) 
Real 


Input:  order  of  MA  model 

input  ;  eoe r f I r 1 ents  of  MA  model 

Input:  variance  of  white  noise 
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UX  Real  Array  (NQ) 

RXO  Real 

IFAULT  Integer 


output:  Rj^(l) . R^(NQ) 

output:  'ly(O) 

output;  failure  indicator  equal  to: 
0  if  no  failure 
1  If  NQ  <  1 


SUBROUTINE  MXCV  (NP,  NQ,  M,  IROWS,  ALPHA,  BETA,  SIGSQ,  RYE,  RYEO,  WKM, 
IP,  RY,  RYO,  IFAULT) 


Formal 

parameters 

NP 

Integer 

Input:  order  of  AR  part  of  model 

NQ 

Integer 

input:  order  of  MA  part  of  model 

M 

Integer 

input:  highest  lag  to  calculate 

(M  >  max  (NP,NQ)) 

IROWS 

Integer 

input:  row  dimension  of  WKM,  IP 

In  calling  program 

(IROWS  >  max  (NP,NQ)) 

ALPHA 

Real  Array 

(NP) 

Input:  coefficients  of  AR  part  of  model 

BETA 

Real  Array 

(NQ) 

input:  coefficients  of  MA  part  of  model 

SIGSQ 

Real 

input:  variance  of  white  noise 

RYE 

Real  Afray 

(NQ) 

output:  Ry^(-I) . Ry^(-NQ) 

RYEO 

Real 

output:  Ry  (0) 

WKM 

Real  Array 

(IROWS, IROWS)  workspace: 

IP 

Integer  Array  (IROWS) 

workspace: 

RY 

Real  Array 

(M) 

output:  R^d),  ....  RyCM) 

RYO 

Real 

output:  Ry(0) 

lFAin.T 

Integer 

output:  failure  indicator 

Failure 

Indications 

Value  of  IFAULT 


Meaning 


0 

1 

2 

3 

4 


no  failure 

NP  -  1  or  NQ  <  1 

M  max  (NP,NQ) 

IROWS  <  max  (NP,NQ) 
singular  matrix  encountered 


The  subroutines  DECOMP  and  SOhV  as  described  by  Moler  (1972)  are  called 
by  subroutine  MXCV. 
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KKSTKICTIONS.  I  IMi;,  STOUACK 

Tf  the  zeros  of  g(z)  are  not  outside  the  unit  circle,  then  one  of 
the  n .  ( j )  will  bo  Rroater  than  or  equal  to  onc>  In  magnitude  thus  giving 
IFAULT  =  10  in  MXPD.  If  the  zeros  of  h(z)  are  not  outside  the  unit  circle 
then  an  element  of  1)^  will  become  nonpositive  thus  giving  IFAULT  =  12  in 
MXPD  or  IFAULT  =  7  In  MAPD. 

The  bulk  of  storage  and  computing  time  in  MXPD  Is  devoted  to  the 
'p  matrix  L^  and  the  x  q  matrix  1.^^.  The  values  and  increase 
as  the  smallest  zeros  of  g(z)  and  h(z)  approach  the  unit  circle. 
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SUUROUr  INE  MXPDl NP , NO. ALPHA.OET A*  SIGSa*Y.IOPT.NritNT2*NHl tNHZ* 
1  NirPO.NPVA«,NPPNH2.  IRON  SI  ,tROMS2.  1  ROM  S  3  .  DEL  .  R  VE  •  R  VE  0.  R  V.  R  VO  •  IP. 
I YWK .RZO.RX.RXO. ALZ I . AL Z . AL X  I • ALX2 .OX • MONE • MT WO . X .E • VPO. PVAR . 

I  IFAU4.T  I 
C 

C  THIS  SUBROUTINE  CALCULATES  PREDICTORS  VPO  AND  lOPTlONALLY) 

C  PREDICTION  VARIANCES  PVAR  FOR  HORIZONS  NHI.....NH2  EACH  FOR 
C  MEMORIES  NTI.....NT2. 

C 

DIMENSION  ALPHAINP  I  .BETAI NO). Y(NT2) .RYEt NO  I  .RVI  IR0WS3). 

1  IP(  I  ROW S3 1  . ALZl I IRQMS3 , IPUWS3) • ALZI IROWSI . NP). Y  WKI NPPNH2I . 

I ALX I ( IRONS3. IROWS3 ) . A L X2 I  I RQWS2 .N 0 1 .DXl IROWS2 I . X( NT2 ) .RX ( NQI . 

1  L( Nr2  I  . VPDI NYPO) .P VARINPVAR ) 

DATA  ZERO. ONE. EPS/ 0.0. 1.0.1 .E-IO/ 

C 

C  CHECK  INPUT  PARAMETERS  : 

C 

I F  A  UL  T  =  1 

IF ( NP.LT. I . OR. NQ.LT. 1 .QR.lOPr.LT. O.OR. lOPT .GT. I  I  60  TO  999 
11  A  UL  T  =  2 

IF  I  NT  1 .LE.NPFNQ.OR .NT2.LT.NTI )  GO  TO  999 
IF  AULT  =3 

IFI NH| .LT. I .0R.NH2.LT.NH1 )  GO  TO  999 
IF  AULT =A 

NCK=(NT2-NT|F1)*(NM2-NH1FII 
IFI NYPO. LT.NCKI  GO  TO  999 
IF  I  NPVAH.LT . 1 »  GO  TO  999 

IFI NPVAR.lt. NCK.AND.lOPT.EQ.l)  CO  TO  999 
IF AULT =5 

IF  I NPPNH2 .LT .NPFNH2 I  GO  TO  999 

IFI IROWSI. LT.NPfZ. AND. lOPT.EQ.I)  GO  TO  999 

IF  I  IROWSI .LT. I  I  GO  TO  999 

IFAULT  =6 

IF  I  IROWSZ.LT.NPFNQFII  CO  TO  999 
IF  AULT  =  Z 

IF  I  IRnwS3.LT .NPfNQ 1  GO  TO  999 
IFAULT  =  B 

IFI SIGSQ.LE.ZEROI  GO  TO  999 
C 

C  FIND  RVO.RVI I  I .. .. .RrlNPFNO)  T 

C 

NPPNO=NPFNQ 

CALL  MXCVI NP.NQ.NPPNQ. I ROWS 3. ALPH A. BE T A . S I CSQ. H YE . R YEO. 

I ALX I . IP.RY.RYO. IFI  I 
IFAULT=9 

IFI IFI. FO. 4)  GO  TO  999 

r 

C  FIND  ALZI.RZO  lALZI  INITIALIZED.  ALPHA13.II  FORMED  IN  ALXllJ.l). 
C  J.I.LE.NP.  RZO  FORMED.  ALZl  FORMED  FROM  ELEMENTS  IN  ALXI)  : 

C 

1FAULT=I0 
DU  20  I=1.NPPNQ 
DO  10  Jsl.NPPNQ 
10  ALZ II  I .3I  =  ZERO 
20  ALZ I  I  I . I  I 30NE 
DO  30  I  -  1  ,  NP 
30  ALXIINP.IMALPHAI  I  I 
C 

IFI NP.FO. I  I  CO  TO  SO 

NPM1=NP-I 

DO  40  J=I.NPMI 


JJzNPMl - J» I 
JJP JJ* I 

PART=AL)(t  (  J  JPI  »  J  JP  I  ) 

IF(PAHT.GE.UNE)  GO  TO  <J9<» 

OEN=:ONE-PART*PART 
no  40  I  si  I JJ 
JJP1MI=J JOI-I 

40  ALX  I  (  JJ  t  I  )  =  (  ALXM  J  JPI  .  I  )-PAflT*Al.X  t4  J  JPI  •  JJPIMI  )  l/OEN 

C 

50  RZ0=SIGSa 

DO  60  Jsl.NP 

60  RZO  =  RZO/(ONE-A|.X1  <  J.J)*ALXI(J«jn 
C 

IF(NP.FQ.I)  GO  TO  80 
DO  70  J=2.NP 
JMl-J-l 

OO  70  Isl«JMI 
J J=JMI-I ♦I 

70  ALZM  J.  IIsALXM  JMI  .  JJI 
80  CONTINUE 
NPP 1 =NPF 1 

OO  90  I=NPPl.NPPNQ 
IFST=I-NPP1 

OO  90  J=1 .NP 
J J= IFST4J 
Jl=NP-JFl 

90  ALZIl  1  .JJ|sALPHA( Jl  I 

C 

C  IF ( lUPT.EQ. n .  find  NONE  ANQ  ALZ  (INVERT  AUZ I »  FIND  ELEMENTS 
C  OF  NP  TM  COLUMN  OF  ALZ  UNTIL  CONVERGENCEt 
C  THEN  FILL  IN  REST  OF  ALZ)  : 

C 

tF( lOPT.EQ.O I  GO  TO  230 
OO  no  l  =  I.NP 
00  100  JslfNP 

too  ALZI I . J)=ZERa 
no  ALZ(I.II=ONE 

IFINP.EQ.t)  GO  TO  140 
OO  130  Js2.NP 
JM1=J-1 

OO  130  Ksl.JMl 
C=ZERO 
JMK= J-K 
JMKPlsJMKFl 

DO  120  IRsJMKPl.J 

120  C=C-ALZ( J.IR)4ALZ(<IR,JNK) 

130  ALZ(JtJMK)=C 
140  CONTINUE 
C 

CK=RZO/SlGSa 
NPP IsNPFI 

ALZINPPI ,NP)3-ALPHA(  1  I 
SUMSQsONEFALPHAI I ) •ALPHAll ) 

MUP=MIN0( IROVS1-NP.NT24NH2-NP) 

OO  160  Js2fMUP 

LLOMsMAXOlO . J-NPIFI 

LUP=  J 

ALOsZEHO 

CC=ONE 

OO  150  LLsLLOW.LUP 


LsLL-1 


no  n  n 


i 


JML - J-L 
L I =NP»L 

IF(L.ar.O)  CC  =  AL2(l.  1  .NP) 

150  ALD=ALr>-ALPHA(  JMLl  *CC 

NPP J=NP» J 
AL?(NPPJ.NP»=ALO 
SliMSQrSUMSQ»ALO*ALO 
IF ( ARSlSU^SQ-CKl.LT.orLl  GU  TO  180 
160  CONTINljF 

IF( MUP.LT.NTZFNHZ-NP)  CO  TO  170 
MUNE=NTaFNH2 
GO  TO  190 
170  CONTINUE 
IFAULT=I 1 
MONE= IRQWSl 
GO  TO  999 
180  MONF^JFNP 
190  CONTINUE 
C 

IF(NP.EQ.t)  GO  TO  220 
NPM 1=NP- I 
Ml MNP=MONe-NP 
no  210  J=t.MlMNP 
NPP JsNP^ J 

OO  210  L  =  I . NPM I 

c=zc«o 

00  200  IR=1.NP 
NPPJMRsNPPJ- IR 

200  C=C>ALPHA( IR )AAL2f NPPJMR.LI 

210  AUZ<NPRJ,L»=C 
220  CONTINUE 
C 

FIND  RXO.RXl 1  I ,•• • tRX (NQI  I 

CALL  MACVlNQ.BETA. SIGSa*RX*RXOt IF  1  ) 

FIND  ALXI.OXI  1  I,... .OXINPFNOI  (NOTE  THAT  ALPHAfl.J)  FOR 
C  I  .LE. J.LE. I.LT.NP  IS  IN  AL Zl ( I » 1 . I- J *  1  I  1  : 

C 

230  continue 

OO  250  I=I,NPPNQ 
no  2*0  J^I.NPPNQ 
2*0  ALXll I .JIs/ERO 
250  ALX 1 ( 1 . I I^ONE 
IFAULT=12 
C 

0X1 1 l=RY0 

0(1  .1*0  (=2*NPPNa 
IMl =1-1 

on  290  J*1 • IMl 
C=ZERO 

00  260  L  =  l  .  I 
00  260  Ms|«J 
ILMM=I ABS(L-M| 

IF(  ILMM.EO.O  >  CsCFALZK  I  «L)  AALZ  I  (  J.M1*RV0 
IF( (LMM.CT.OI  CsC*ALZIIl*L}#ALZI< JvMIARYl ILMN) 
260  CONTINUE 

IF (J.EO. 1 )  GO  TO  280 
JMl *J-l 

OO  270  L«I.JM1 

C=C-*LXI  I  1  .L  l*OX(L)*ALXl(  J.L) 


270 


?80  ALXI ( I . J ) =C/OX( J) 

290  CONTINUE 

C=RX0 

IFII.GT.NPI  GO  TO  320 
C=ZERO 

OO  310  . I 

CI  =  A1.ZI  (  I  .Ml 
C2-ZeM0 

00  300  L=l.I 
ILMMzI ABS(M-L) 

IF(  ILMM.EQ.OI  C23C2*’ALZt<  I  .LI  AM  TO 
IFIILMM.GT.O)  C2=C2*ALZI( 1*  LI*RV(  ILMM I 
300  CONTINUE 

310  C=CFCIAC2 

320  DO  330  L=t*IMl 

330  C=C-OX(Ll  AALXl  (  1 ,1.  I  «  ALX  t  <  t  .L I 

IFIC.LT.EPS)  GO  TO  999 
340  DXIII^C 
C 

C  FIND  MT«a.  ALX2.  AND  THE  REST  OF  OX  (IN  THIS  SECTION* 

C  I  AND  . . t-l  REPRESENT  THE  INDICES  OF  THE  NONZERO. 

C  NONONE  ELEMENTS  OF  THE  MATRIX  LSUBX.  NOTE  THAT  |F  INTEGERS 
C  M.LE.N.LE.NPFNQ*  THEN  LSUBXIN.M}  IS  STORED  IN  LSUBX(N.M) 

C  WHILE  IF  N.GT.NPANQ.  then  LSUBXCN.MI  IS  IN  AL X2(  N-NP-NO . 

C  NO-(N-M)^I).  M=N-NQ« ... ,N- II  : 

C 

IFAULT=I2 

IUP=MIN0( IROWS2-NP-NQ. NT2ANH2-NP-N0I 
OO  400  list .IMP 
l=NPPNQAl I 
IMNO=|-NQ 

ALX2(  I  I . 1  I^RXI NQI/OXI  IMNQl 
IF(NQ.EO.l)  GO  TO  370 
OO  360  JJsZ.NQ 
J= IMNG* J J-1 
NOI^NQ- JJFt 
Jl = J-NPPNQ 
C=RX( NQI ) 

JJMI=JJ-1 

OO  3S0  LL^t.JJMI 
L - IMNOFLL-I 
J2=NQ-( J-LIFI 

IF!  J.LE.NRPNQ)  CtsALXKJ.LI 
IF( J.GT, NPPNO)  C I=ALX2(3t *32) 

350  C=C-ALX2< II.LL)40X<Ll4CI 

360  ALX2(  I  I . JJ|3C/0X( J  I 

370  C=RXO 

DO  3B0  L-l .NO 
LL* I-NQAL-I 

380  C-C-ALX2(  II  .L)*ALX2(  I  I.LMOXILLI 

OX(|)=C 

4F( OXI I I.LT.EPSI  GO  10  999 
C 

DO  390  JJ^l.NQ 
NOMJPIsNO-J J4| 

390  IFIABS(ALX2(  II.3JI-BETA(NQMJP|  M.aT.OEU)  GO  TO  400 

IFI ARSIDXI I l-SIGS0).GT.0ELI  GO  TO  400 
MTWO= I  I ANPPNQ 
GO  TO  420 
400  CONTINUE 

IFI  IUP.LT.NT2ANH2-NPPNQI  GO  TO  410 


n  n 


Mrwa=NT2»NH2 
GO  10  420 
0  CONTINUE 
IF  AULT=  I  J 
MTITO=  IHOWSZ-NPPNQ 
GO  TO  0<J9 

>  CONTINUE 

FIND  X<1) . X( NT21 .El  1  )  .Ef Nrai 

NPP I =NP*I 
X<  I  )  =  V<  I  I 

IFCNP.EQ.I)  GO  TO  450 
on  440  J=2«NP 
C-Y<  Jl 
JMl-J-l 

DO  430  L=1.JMI 
JML  = J-L 

)  C=CFAUZI  (J.  JMLI  •VI  JML) 

)  X|J)=C 

)  OU  470  J=NPPl.Nr2 
C=V( J) 

DO  460  L=1.NP 
JML= J-L 

)  C=C4ALPHAIL )*Vl JML  ) 

I  X|J)=C 

El  1  |3X(  1) 

OO  490  Ja2.NPPNQ 
C=X ( J ) 

JMl =J-I 

00  480  Lsl.JMI 
)  C*C-ALXII J.L)*elL) 

»  EIJ)=C 

MLOW  =  NPnNQF  t 
MUP^sMINOINTZ.MTMO) 

OO  510  JsMLOM.MUP 
JJ= J-NPPNO 
C=X I J ) 

00  500  L=l .NO 
LL=J-NQ»L-l 

I  C=C-ALX2I JJ.LIAEILLI 

>  EIJ)=C 

IFI MUP.EQ.NTZ)  GO  TO  540 
MUPP1=MUPF1 
OO  530  JsMUPPt.NTa 
C=Xl J) 

00  520  LSI .NO 
JMLsJ-L 

CsC-BETAILMEI  JMLI 

•  E(J)sC 

>  CONTINUE 


INO  vpo  : 

NP0PTSNH2-NMI4I 
00  630  NTsNri,NT2 
NSUFAM=INr-NTl )*NPOPT 
NTMNPsNT-NP 

OO  550  |s|.NP 
I  I sNTMNPAI 


550 


T 

-  ^ 

VHK(  I) =Y(  I  I  I 
DO  510  NH=1,NM2 
NPPNH=NP^NH 
NTPNH  =  NT4-NH 
IROV»LX  =  NTPNM-NPPNQ 
XTPH=ZERO 

IF(NH.aT.NO)  GO  TO  590 
IF( NTPNH.Gr«MTWO)  GO  TO  570 
00  560  K=NH.NO 
INOL=NQ-K* I 
INDE=NTPNH-K 

XrPH=XTPHf AL X?( 1R0MLX«1N0L) *E( INDE) 

GO  TO  590 

00  580  K=NM.NQ 
lNOE=NrPNH-K 

XrPH=xrPM*-aETA(K)*E(  JNOE) 

C*XTPH 

DO  600  J=l.NP 
INOY=:NPPNH-J 
C=C-ALPHA( J I *YMK< INDY  1 
YVIKCMPPNH)  =C 
DO  620  l=NHI.NH2 
N1 0=NP* I 
N1  t  =  I-NHt  F  t 
N1 2=NS0F ARFNI 1 
YPO(N12) 2YWK(NI0) 

630  CONTINUE 
C 

C  IF  I0PT.E0.1«  PINO  PVAR  : 

c 

IF( lOPT.EO.O)  GO  TO  690 
Ml MNP=MONE-NP 
DO  680  NT=NTI.NT2 
NSOFARsINT-NTl l^NPDPT 
DO  680  NM=NHI«NH2 
NINOX=NSOFARFNH-NHI ♦! 

NTPNH=NTFNH 
C=ZERO 

DO  670  KPl=lfNH 
K=KPI- 1 
KROW=NTPNH-K 
C1=SIGSQ 

IFf KROW.LT.MTMO)  CIsOXfKROWl 
C2=aNE 

IFIK.EO.O)  go  to  670 
NPPK=NP»K 
CZaALZINPPK, NPI 
DO  660  |R«I.K 
KNIR  =  K-IR 
C3=ZEPa 

IF ( KMR.GT .MIMNPI  GO  TO  640 
INDL=NPfKMR 
C3=Al.Z<  INOC.NPI 
640  C4=7EPn 

IFfIR.GT.NO)  GO  TO  660 
N10=NIPNH-KMR 
IFINIO.GT.MTNOI  CO  TO  650 
N1 I =Nl O-NPPNO 
N1 2=NO-IHFI 
C43ALX2(N1 I.N12) 

GO  TO  660 


560 

5  70 

580 

590 

600 

610 

620 


6bO  C*  =  fJETA(IH) 

ftftO  C?-C2*C3*CA 

670  t=C*Cl*C2*C2 

680  PVAW(NINUX ) =C 
C 

690  continue 

c 

ir  AUL  r^o 
999  MET  URN 
END 

SUBROUTINE  MAPOINQ .BETA. SIGSO. V. lOPT.Nri.NTa.NHl.NHZt 
INYPO.NPVAR*  IHOWSt 0EL«RX.RX0  rOX • AL X, HT WO .E • YPO* P VAR • 1 F AULT I 

C 


C  THIS  SUBROUTINE  CALCULATES  PREDICTORS  YPO  AND 
C  (OPTIONALLY!  PREDICTION  VARIANCES  PVAR  FOR  HORIZINS 
C  NHI.....NH2  EACH  FOR  MEMORIES  Nri*...*Nr2. 


C 

DIMENSION  BErA(NQ) « Y(Nr2 ) .RXfNQI. DXI I ROWS I . AL X( 1 RO WS. NO ) • 
IF( NT2 I . YPD( NYPO> .PVAR  <  NPVAR ) 

DATA  ZERO .ONE .EPS/0.0. I. 0. t .E- I 0/ 

C 

IFAULT=I 

IF( nQ.LT . 1 .OR. iOPT .LT. O.OR. lOPT.cr.l )  GO  TO  199 
IFAULT-sZ 

IF  (  NTI  .lT.NQ*-!  .OR.NTl  .GT.NT2I  GO  TO  199 
IF AULT=3 

IF( NH| .LT . I .OR.NHI .GT.NH2)  GO  TO  199 
1FAULT=4 

NCK  =  <NT2-Nri  PI  lAINHZ-NHlf'l  I 
IFINYPO.LT.NCKl  GO  TO  199 
IFINPVAR.LT. I )  GO  TO  199 

IFiNPVAR.LT.NCK. AND. IOPT .EQ.I)  GO  TO  199 
IFAULT=5 

IF ( IRUWS.lt .21  GO  TO  199 
IF  AULT  =  6 

1F< SIGSQ.LE.ZEROI  GO  TO  199 
IFAULT=7 
C 

C  FIND  MTWO  AND  ROWS  OF  AL X .  OX  : 

C 

CALL  MACVINO.BETA. SIGSQ.RX.RXO. IF  I  I 

C 

DXl 1 )=RX0 
DO  10  1=1. NO 
ALX( I. I |=ONE 
E(  1  1=Y(  I  I 

lUPsMINOI IR0WS.NT2  ANH2 I 
on  100  132. I UP 
I  IsMAXOl  I-NO-1  .01*’  I 
IMI3I-I 

NELTS= IMl- I  1*1 

DO  30  J3|1.IM1 
JINOaJ-l 1*1 

IMj3l - J 

JMlsJ-l 
C=RX(  IMJ  I 

IFiJ.EO.Itl  GO  TO  30 
Jl 3MAX0( J-NO-I .01*1 
DO  20  LsII.JMI 
LL=L-I 1*1 
JJ»L- Jl* I 


10 


r>  n 


20  C=C-ALX  I  I  ,LL  I  *f>X  <LI*ALX<  J.  J  J  I 

30  ALX» l.JINDI-C/DX< J» 

C=RX0 

Ua  40  J=II.IM1 
JIN0=J-I I *l 

40  C=C-OX( J )*ALX< I • JIN0)*ALX( ( . Jf  NO) 

IF(C.LE*EPS)  GO  TO  199 
DX< I )=C 


IF( 

i.gt. 

NT2  ) 

GO 

TO 

60 

C=Y 

(  I  I 

DO  50 

J=  1 

•NELTS 

11  =  1- 

NFLT 

SF  J 

-1 

50 

C=C-ALX(  1 

.  J) 

•El 

I  I  ) 

F(  I 

)  =C 

IF( 

I.LE. 

NO) 

GO 

TO 

1  00 

60 

DO  70 

J=1 

•  NO 

• 

J  J=NQ 

-  JF  1 

70 

IF  1  ABS(  AL 

X(  I 

.  J  ) 

-BETA! JJ ) l.GE.OEL 1 

IF( ABSIOXl I )-SIGSQ l.GE.OEL)  GO  TO  100 
MTWn= I 

IF(I.GE.NT2)  GO  TO  110 
IPl =1 ♦ 1 

DO  90  J=IP1.NT2 
C=Y( Jl 

DO  80  K*1.NQ 
JMK=J-K 

80  C=C-BETA|K)4C< JMK) 

90  E1J)=C 

GO  TO  120 
100  CQNTINVIE 
IFAUUT=0 

IF<  IUP.I-r.Nr2*NH2l  GO  TO  199 
110  MTWO=NT2FNH2 
120  IFAULT  =  0 

CALCULATE  PREDICTORS  : 

NPOPT=NH2-NHl ♦ 1 
DO  140  NTsNTl ,Nr2 
NSOFAR=lNT-NTl )4NP0PT 
DO  140  NH=NHl<NH2 
NIN0=NS0FAHFNH-NH1  1 
YPO<NINO)=ZERO 
IFINH.GT.NQ)  go  to  140 
NTPNrt=NT4NH 
C=ZERO 

DO  130  K=NHtNQ 
I NOE=NTPNH-K 
INI)L  =  NO-K4l 
C1=6ETA1K» 

IF( NTPNH.LE. MTWO)  C 1 = ALX ( NT PNH» I NOL » 
C  =  CFC14E1  INUEI 
YPD1NIN0)=C 

IF  lOPTsl.  CALCULATE  VARIANCES  : 

IFI lOPT .EO. 0)  GO  TO  199 
DO  180  NT=NTI,Nr2 
NSOF AR=INT-NT1 )4NPOPT 
DO  100  NHsNHI*NH2 
NHMIsNM-I 


130 

140 

C 

C 

c 


too 


n  n 


NINt)-NSOFAR*NH-NHJ  ♦  1 
NTPNM-Nr ♦NH 

c=s  ir.so 

IF  (  NT  .L  T  .MTWin  C^OXTNTPNHI 
IF ( NH.EQ . 1 )  GO  TO  I flO 
NHUP-MINOINHMl .NO) 

IFCNT.GE.MTWOI  GO  TQ  160 
on  150  K-l.NHUP 
I  Nr>D»NTPNH-K 
IN0L=N0-KF1 

t  50  C=:C*OX<  INOO  I  ♦ALK  TNTPNH.  INOL  »*ALX(  NTPNH,  I  NOLI 

GO  Tu  lao 

160  on  1/0  K=1,NHUP 

1/0  c-CFS iGSQ*BerA(K )*Ber ACK> 

IflO  pvar(ninoi=c 

c 

199  PFTURN 
FNO 

^UnROUT INF  ARPOl NP. AL PH A . S I GSO • V * I0Pr.NTI.NT2*NHl cNHa* 

I NYPn.NPPNH2 . YMK .GAM«yPO. PVAR. IFAULT ) 

C 

C  THIG  SUHROUTING  CALCULATES  PREDICTORS  YPO  ANO  lOPTlONALLrl 
C  PREDICTIUN  VARIANCES  PVAR  FOR  HORIZONS  NHU....NH2  EACH 
C  FOR  MEMORIES  Nri.....Nr2 

c 

DIMENSION  Y<NT2I.ALPHA(NP| . VWK I NP PNH2 I . GAM C NPP NH2 ) • 
lYPOI NYPOl .PVARINH? I 
DATA  ZERO.ONF/0.0. I .0/ 

C 

IFAULT=1 

IFlNP.LT.n  GO  TO  100 
IFAULT  =2 

IF ( NTl .LE.NP.OR.NT I .GT.NTZI  GO  TO  100 
IF AULT=J 

IF  I NHl  .LT. 1 .OR.NHl  .GT .NH2I  GO  TO  100 

IF  AULT-4 

1F( SIGSQ.LE.ZERU)  GO  TO  100 
IF AULT=5 

ir( lOPT.LT^O.OR.iaPT.GT. I )  GO  TO  100 
IFAULT=6 

NCK^I NT2-NT 1 ♦! I •(NH2-NHI ♦! I 

IF! NYPO.Lr.NCK.OR.NPPNH2«LT.NPFNH2l  GO  TO  100 
IFAULT=0 
C 

FIND  PREDICTIONS  : 

NPf)PT  =  NH2-NHl  1 
no  50  NT=NT1.NT2 
NSOF AR=(NT-NT I »*NPOPT 

ntmnp=nt-np 

DO  10  1=1. NP 
I i=Nr-NP* I 
10  VMKI  I  )  =  Y(  I  I  I 

DO  30  NH=1,NH2 

nppnm=npfnh 
ntpnh^ntenh 

C=:ZERO 

on  20  I-»l  .NP 
I I=NPPNH-I 

20  C=C-ALPHA(  I  I  •YIIKII  I  » 

30  YWKINPPNHisC 


DO  40  NH-NH1.NH2 
INONH  =  NSOF4R«-NH-NHI  ♦! 

INOV*K=NP«-NH 

40  YP0(  INONHI-sVWKl  INDMKI 

50  CONTINUE 

IF  iaPT=l.  FIND  VARIANCES  : 

1F( lOPT.EQ.O)  GO  TO  100 
GAM(  1  ) =-ALPHA(  1  I 
|F(NH2.Ea.l)  GO  TO  80 
DO  70  NH=2.NH? 

LLOW=MAxa(0,NH-NPl  «- 1 
C=:2ERO 

00  60  LL  =  L1.0W.NH 
L=LL-1 
C I^ONE 

IFIL.GT.O)  CI=GAM(LI 
NHML=NM-L 

60  C=C-ALPHA ( NHML » *C1 

70  GAN(NH)=C 
80  PVARl  1  l-SIGSO 

IF(NH2.Ea.ll  GO  TO  100 

00  90  l=2«NH2 

IM|=I-1 

90  PVARI I l^PVARI |M| |4SIGS04GAN( IMl l*GAN( I Ml  1 
100  RETURN 
END 

subroutine  MAtV<NQ«0ETA«SIGSQ.Ry»RyO*lFAULT» 

THIS  SUBROUTINE  CALCULATES  MA4NQ)  AUTOCOVAMI ANCES  OF  LAGS 
0....«Na  (NQ.GT.OI 

DIMENSION  BE TAl NQI .RTl NO) 

DATA  ONE/I .0/ 

IFAULT=1 

IFINQ.LT.ll  GO  TO  40 
IFAULT*0 

C=ONE 

OO  10  l-l«NQ 
10  C=C*BEr A( 1 lABETAI  I  ) 

RyO=C4SIGSQ 

DO  30  IVsl.NO 
C=BETA( IV> 

ir( IV.EO.NQ)  GO  TO  30 
NQMI  V=N(1-IV 

OO  20  J^I.NQMIV 
JPI V= J* I  V 

20  C-C»8ETA<  JMBETAl  JPIV) 

30  RTl 1V)=C4$IGS0 
40  RETURN 
END 

SUBROUTINE  MXCVf  NP  .NO*  M,  IROliS*  ALPHA*  BE  T  A.  S I GSQ  tliyE  (RyEO  • 
I «KM  * IP.Ry.RVO. IF AULT I 


THIS  SUBROUTINE  CALCULATES  ARMA(NP*NQI  AUTOCOVARIANCES  FOR 
LAGS  0.****M  1 M.GE.MAXl NP.NQI •  NP*NQ*GT*OI  i 


OIMnNSIDN  ALPHA<NP) . BET A f NQ I *9 YE( NQ I t WKMI IRONS. IRONS)  • 
1 IP(  IRONS) .R Yf  M) 

DATA  ONF.2ERO/ 1 .0. 0.0/ 

C 

IFAULT=  1 

IF< NU.LT.I .OR.NP.LT.l )  GO  TO  110 
IFAULT=2 

MAXPO  =  NIAXO(  NP.NQ) 

MM3HAXPQF1 

IF( M.LE.MAXPQ)  GO  TO  110 
IFAULT=3 

IF ( IRONS. Cr. MM)  GO  TO  110 
IFAULT=A 

C 

C  FIND  RYEO.RYEl  1  I . . • • .RYE <NQ)  : 

C 

HYF0=SI GSO 
OO  30  IVxl.NQ 
C  =  S IGSOABETAl I  V) 

NUP=M INOI I  V ,NP| 

DO  20  J= 1 .NUP 
1 VMJ-I V- J 

IFl tVMJ.EQ.O )  GO  TO  10 
C =C -ALPHA! JJ*RY El IVMJ  ) 

GO  TO  20 

10  C=C-ALPHA( JlARYEO 

20  CONTINUE 

30  RYe<IV)=C 
c 

C  USE  OECOHP.  SOLV  TO  OBTAIN  R YO.RY f 1 )••••• RY! MAXI NP • NQ > ) 

C 

00  40  |V=1.MM 
RY I  IV  )  =  ZERO 

OO  40  J-l.MM 
40  NKMI IV. JlsZERO 

C 

NPP  I  =  NP«- 1 

NOP 1=N0F I 

no  60  IVPI =1 .NQPl 

IV=IVP1-1 

C=RYF0 

irilV.GT.O)  C=C48ETA|IV) 

IFI  IV.EQ.NQ)  GO  TO  60 
OO  SO  K^IVPl.NQ 
KMI V=K-I V 

SO  C=C4HeTAIK)4RYelKMI VI 

60  RYI I VPl )=C 
C 

DU  70  IVP1=1.MM 
I  Vs IVPl-1 

NKMI  I VPl . I VPl  )=NKMI  I  VPl. IVPl |♦ONE 
DO  70  Jsl.NP 
I  1^1 ABSl  IV-JI41 

70  NKMI IVPl . t 1 )sWKM| 1 VPl. I 1)4ALPHAI J ) 

C 

CALL  OECUMPl MM. IRONS. NKM.IPI 
IFI  IPIMMI.EQ.Ol  GO  TO  no 
IFAULTsO 

CALL  SOLVIMM. IRONS. NKM.RV.IPI 

RYO=RY( I ) 

on  80  IVsl.MAXPO 


lVPI=IVf 1 

80  RV<  IV)=R>r(  IVPl  ) 

c 

C  USE  DIFFERENCE  EQUATION  TO  GET  THE  REST  OF  THE  RT  : 
C 

IF( M.EQ^MAXPQ)  60  TO  110 

DO  100  IV=MM.M 

C=ZERO 

DO  90  Jsl.NP 
I VM J=IV- J 

90  C»C-ALPHA( J)*RYf IVMJ) 

too  RYC  IV  »  =  C 
C 


110  RETURN 
END 


END 
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